Quantum Simulations of Relativistic Quantum 
Physics in Circuit QED 



J S Pedernales 1 , R Di Candia 1 , D Ballester 21 and E Solano 1,3 

1 Departamcnto de Quimica Fi'sica, Universidad del Pais Vasco UPV/EHU, Apartado 
644, 48080 Bilbao, Spain 

2 Department of Physics and Astronomy, University College London, Gower Street, 
London WC1E 6BT, United Kingdom 

3 IKERBASQUE, Basque Foundation for Science, Alameda Urquijo 36, 48011 Bilbao, 
Spain 

^ Abstract. We present a scheme for simulating relativistic quantum physics in circuit 

quantum electrodynamics. By using three classical microwave drives, we show that 
a superconducting qubit strongly-coupled to a resonator field mode can be used to 
simulate the dynamics of the Dirac equation and Klein paradox in all regimes. Using 
the same setup we also propose the implementation of the Foldy-Wouthuysen canonical 
transformation, after which the time derivative of the position operator becomes a 
constant of the motion. 



Quantum Simulations of Relativistic Quantum Physics in Circuit QED 



2 



1. Introduction 

R. P. Feynman is credited for introducing the idea of a quantum simulator [1]. At 
his keynote speech during the 1st Conference on Physics and Computers in 1981, he 
claimed that because the laws of Nature are not those of classical mechanics, it would 
be better to build a quantum simulation of a physical problem under the same quantum 
mechanical laws. For the following ten years, his contribution was somewhat hindered 
by the work of Deutsch [2] , who showed the existence of a universal quantum computer, 
a quantum mechanical generalization of the classical Turing machine using quantum 
logic gates and algorithms. 

The concept of quantum simulators will reappear later in 1996, when Lloyd 
presented a general proof for Feynman's conjecture [3]. Quantum simulators can surpass 
the capabilities of classical ones, where the complexity of a problem scales exponentially. 
Lloyd argues [3] that whereas the simulation of a general 40 spin— 1/2 particle system 
would be enough to outperform existing classical computers, the factorization of a 100- 
digit number with Shor's algorithm would need of thousands of qubits to become 
nontrivial. Therefore, it is likely that the first quantum device to beat classical 
computers would be a quantum simulator. But there is another reason why quantum 
simulators are attracting so much attention lately. They provide powerful analogies 
to understand the physics of most varied type of systems one can think of, including 
condensed-matter physics, high-energy physics, cosmology and quantum chemistry [HE]. 
When compared to the original systems, quantum simulators show an unprecedented 
degree of controllability over all physical parameter regimes. 

In this work, we show how the physics of relativistic quantum mechanics can be 
simulated in a very different setting. That of circuit quantum electrodynamics (cQED) 
and superconducting qubits. We will focus our attention on the fundamental equation 
in relativistic quantum physics, the Dirac equation. Proposed by Dirac himself in an 
attempt to combine quantum mechanics and special relativity in 1928 [6], this equation 
is the first one to describe elementary spin— 1/2 particles — such as the electron — and 
the existence of antimatter. Here, we restrict our discussion to the simplest case of 1+1 
dimensions for which the algebra of the Dirac spinors can be reduced to that of standard 
Pauli matrices [7j. In the standard representation [7], the 1+1 Dirac equation can be 
written as 



where a y and a z are Pauli matrices, P represents the momentum of a Dirac particle of 
mass m and c the speed of light. Soon after its introduction, Schrodinger realised that 
there is something odd in the solution of this equation [8]. The position for a free Dirac 
particle can be integrated to give [9] 




(1) 



x{t) = x + c 2 PH~H + Z(t) 



(2) 



Quantum Simulations of Relativistic Quantum Physics in Circuit QED 



3 



with H = mc 2 a z + cPa y and 

Z(t) = l -hc (a y - cPTT 1 ) H- 1 (exp (-i2Ut/h) - 1) . (3) 

This third term appears because the time derivative of the position operator is not 
a constant of motion, i.e. [x,"H] 7^ 0, even for a free particle. Given a general 
initial wavepacket, the expectation value (Z(t)) shows an oscillating pattern with 
amplitude ~ h/mc and frequency ~ 2mc 2 /h that is known as Zitterbewegung or jittering 
motion [9]. This behaviour stems from the interference between positive and negative 
energy components. Indeed, it is well established that initial states with only positive or 
negative components will not show any Zitterbewegung [ZJ E]- In spite of the appealing 
physics behind the Dirac equation, there has never been an experiment to test the 
existence of such an effect. In fact, for relativistic electrons, one can estimate a 
Zitterbewegung with an amplitude of about 10 -4 nm and a frequency of 10 21 Hz, making 
its eventual detection very demanding. 

There is another interesting property of the Dirac equation related to the behaviour 
of relativistic Dirac particles under the effect of an external scalar potential, 

inS- = {mc 2 a z + cPa v + <$>{X)}il). (4) 

In case of a nonrelativistic particle impinging against a semi-infinite step potential $(A), 
we know that if the kinetic energy of the particle is smaller than the potential height, 
no propagating solutions exist at the other side of the potential barrier. In other words, 
the particle penetrates slightly, but it will be eventually reflected. Klein found out 
that Eq. (jlj) for a relativistic Dirac particle allows for propagating solutions beyond the 
potential barrier [TO] . And in the case of an infinitely high potential, the tunneling occurs 
with unit probability. More precisely, while for a small potential $ < E—mc 2 the particle 
with kinetic energy E will be mostly transmitted, for a larger one E—mc 2 < $ < E+mc 2 
it will be reflected. However, making even larger the potential $ > E + mc 2 the particle 
will propagate through the barrier. To explain this apparent paradox, one can make 
use of the idea of particle-antiparticle creation by the potential [7]. A particle can go 
through the potential by turning into its antiparticle. 

Although the Dirac equation is considered a milestone of relativistic quantum 
physics, playing a crucial role towards the later development of quantum field 
theories, in the last decades it has driven a renewed interest in different playgrounds. 
Examples include: superconductivity, where Bogoliubov quasi-particles can exhibit 
Zitterbewegung [TTj : semiconductor physics, with a close analogy between the k ■ p 
theory of energy bands in narrow-gap semiconductors and the Dirac equation [12]; the 
physics of graphene, characterised by quasi-electrons behaving like relativistic Dirac 
particles in two dimensions [13J; a proposal for generating relativistic scattering in optical 
latices [29J ; photonic arrays with implementation of the quantum Rabi model and the 
Dirac equation [301 IS]; and the physics of ion traps, where the vibrational and internal 
degrees of freedom of the ions can be made to follow relativistic quantum dynamics 
[T4l fT5l fT6l [T71 [18] . or even behave as interacting bosons and fermions j32j [33] . In the 
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spirit of quantum simulations, the latter represents an important advancement in terms 
of full quantum controllability of the dynamics and the retrieval of information through 
measurement. 

Here, we present a scheme for performing quantum simulations of relativistic 
quantum mechanics using state-of-the-art cQED technology This is motivated by 
two reasons. First, since 2004 [T9| [20] . cQED has become the quantum platform 
with the most promising perspectives in terms of scalability and coherence. Second, 
there are crucial physical differences with respect to any previous implementation of the 
Dirac equation and Klein paradox. In our case, the physics of a relativistic spin— 1/2 
particle is simulated by the interaction of two degrees of freedom of different physical 
systems, i.e., a standing wave in a superconducting resonator and the two lowest 
levels of a superconducting qubit, where none of them move at all. The position and 
momentum of the Dirac particle are then encoded in its phase-space or field quadrature 
representation. This opens up new possibilities for combining intracavity fields with 
propagating quantum microwaves [211 E2] in scalabale quantum network architectures 
with delocalised and/or sequential interactions. 

2. General derivation 

For our protocol, we require a two-level superconducting qubit, working at its symmetry 
point, strongly coupled to a single electromagnetic field mode of the resonator. This 
interaction will be described by the Jaynes-Cummings model (JCM) [231, [191 I2Q]. In 
addition, we introduce three external classical microwave drivings, two transversal to 
the resonator [21], coupling only to the qubit, and one longitudinal that only sees the 
resonator mode. The time dependent Hamiltonian of the complete system is given by 

U = ^cr z + hwtfa - kg (Va + atf) - HQ (e^ t+ ^a + e^^M) 

- h\ {e i{ut+v) a + e- l{ut+ ^o- ] ) + ^ (e^a + e^tf) . (5) 

with a y = ia — ia^ = i\g)(e\ — i\e)(g\ and a z = \e)(e\ — \g)(g\, \g) and \e) being 
the ground and excited states of the qubit. Here, Tko and hui q stand for the photon 
energy and qubit energy splitting, whereas g is the coupling constant. Likewise, the two 
orthogonal drivings have real amplitude Q, A, phase ip, and frequencies u, v. Finally, 
the longitudinal driving is characterised by an amplitude £ and a frequency u. Note 
that while we have chosen two of the drivings to be resonant with the resonator field 
mode, the other parameters will be set later on. For simplicity in our presentation, we 
will also assume that u q = u, i.e. qubit and resonator field interact resonantly as well. 

Our derivation consists of two straightforward transformations. First, Hamiltonian 
in Eq. (jHJ) can be simplified using the reference frame rotating with the resonator 
frequency u 

U Ll = -Hg (a^a + aa f ) - hVl (e itp a + e^V) + hi (a + a f ) 

- h\ (e i K"- w ^ a + e -<[("-«)^] o-t) . (6) 
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Second, we will write this Hamiltonian in another interaction picture with respect 
to the time independent term T-Lq 1 = —Ml [e %ip o + e~ l<p o~^\ . The resulting expression 
reads 

^=-y ({!+)(+! -\-)H + e- l2Qt \+)(-\ 

- e i2m |-)(+|} e^a + H.c.) 

- y ({l+X+l - I — >< — I - e ' l2m l+)H 

+e i2m |_)( + |} e i(—)t + H.c.) + fif (a + a f ) , (7) 

where we have introduced the rotated spin basis |±) = (|#) ± e~ iv \e)) /y/2. To simplify 
this Hamiltonian expression further, we can now set u — v = 2Q, and assume the 
amplitude of the first driving Q to be large as compared to the rest of frequencies in ([7]). 
This implies we can apply the rotating-wave approximation (RWA) to yield the effective 
Hamiltonian 

h\ kg „ . ^ r- A 
rt cfr = — cr 2 + -yi "^ + 2 x, (8) 

where ip = n/2 and we make use of standard quadratures of the electromagnetic field, 
i.e. dimensionless x = (a + a)) /y/2, p = —i(a — al)/\/2, at variance with X and P in 
Eq. (TJJ, with the commutation relation [x,p] = i. 

The Schrodinger equation of our system resembles that of the 1+1 Dirac equation as 
in Eqs. (EQ) and 01]), where the terms hg/\/2 and h\/2 are related to the speed of light and 
the mass, respectively. In addition, we have an external linear potential $ = h£\/2x, 
that depends linearly on the position of the particle. Clearly, in the simulated dynamics, 
one can cover a wide range of physical parameters. While, for a fixed coupling strength 
g, the simulated mass is proportional to the amplitude of the weak orthogonal driving 
A, the strength of the linear potential can be tuned through the amplitude £ of the 
longitudinal driving. This is an interesting advantage over the implementation in ion 
traps, where a second ion is needed to simulate the external potential [16j \T7\ . 

As already mentioned, the study of relativistic quantum effects, such as 
Zitterbewegung or Klein tunneling, should be done in the phase-space representation 
of the electromagnetic field in the superconducting resonator. In all numerical cases 
studied below, the initial state of the bosonic degree of freedom of the Dirac particle is 
assumed to be a wavepacket with average position (x ) and average momentum (p ), 

iP( x ) = 7T- 1 ^exp{t(p )x} exp /_ ~ fr)) 2 1 (Q) 



2 

This coincides with the x— quadrature representation of a coherent state of the 



electromagnetic field 



(£oH|£o>^ = x> ( >°>+£W ) |0), |0) being the vacuum state of the 

bosonic field, and T>[a) = exp {aa) — a*a} the coherent displacement operator. 

For the numerical simulations discussed below, we have chosen a value of the 
superconducting qubit and resonator frequencies of u q = uj = 2n x 9 GHz, a qubit-field 
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coupling of g = 2tx x 10 MHz and a strong transverse driving amplitude of f2 = 2tx x 200 
MHz. The interaction time needed is of about 60 nsec, which is well below standard 
decoherence times in cQED experiments. 



3. Free Dirac particle 

In the absence of external potential, $ = 0, the Dirac equation for a free particle ([I]) 
does not couple the different spinor components. The Dirac Hamiltonian 

H D = —a z + -j=£yP (10) 

has two important limits depending on the value of the mass of the particle. First, in the 
massless case the Hamiltonian reduces to Hd = hg/y/2a y p allowing for straightforward 
interpretation of its dynamics in terms of phase-space representation. Fig. Ufa) shows 
the evolution of a massless particle with initial state |+, 0). Being its spin in an eigenstate 
of a y , the effect that the interaction has is just to generate a coherent field displacement 
from its original position, the vacuum state |0), to the final one in which the state of 
the field is \gt/2). This simulates the evolution of a massless free Dirac particle moving 
at a speed of light stemming from hg/y/2. As we expect, the dynamics of this massless 
particle would be essentially unchanged if one assumes the initial state to have a non- 
zero value of the expectation value of momentum. This is depicted in Fig. QJd), where 
the particle starts its trajectory with (po) = 2. The particle covers exactly the same 
amount of space in the x— quadrature, ending in the coherent state \gt/2 + y/2i). 

Second, when the mass of the particle becomes large enough, the dynamics should 
tend towards that of a non-relativistic one evolving under a free Schrodinger equation. 
For h\/2 > hg/V2(p), i.e. mc 2 > c(P), one can derive a second-order effective 
Hamiltonian [H] for Eq. flS} to yield H^Rei = o~ z H(g 2 / X)p 2 . This is indeed the 
nonrelativistic Schrodinger Hamiltonian we expected, but with the addition of the spin 
operator o~ z , which is reminiscent of the spinor- momentum coupling of the original Dirac 
equation. This spin operator ensures the existence of two decoupled equations, valid for 
positive and negative kinetic energies. In Fig. Hfc), we depict the evolution of the initial 
state |e, 0) under the time dependent Hamiltonian ([6]) for the case in which the mass of 
the particle is such that hX/2 = 4 x hg/y/2. This generates a squeezed vacuum state 
with degree of squeezing increasing linearly in time. This behaviour is in agreement with 
the aforementioned nonrelativistic approximation for the Hamiltonian, proportional to 
p 2 . Such a simulation realizes a standard example in all quantum mechanics textbooks. 
Namely, the free evolution under the nonrelativistic Schrodinger equation of an initial 
wavepacket with (po) = produces a particle that remains centered at the same initial 
position, with a wavepacket spreading over time in the x— quadrature. However, the 
effect of the operator p 2 is not simply squeezing. This can be clearly seen in Fig. [Tff) 
where we have taken a particle with initial state |e, \/2i), that is a wavepacket with 
non-zero kinetic energy such that (po) = 2. Now the wavepacket not only spreads over 
time, but its centre moves linearly in time to the right as the expectation value of the 
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Figure 1. Wigner function W(x,p) representation of the state of the electromagnetic 
field mode inside the resonator. The initial state evolves under the Hamiltonian © 
for a time period of 60 nsec. The physical parameters used are g = 2tt x 10 MHz, 
fl = 2tt x 200 MHz, f = Q, together with: (a) A = being the initial state |+,0>; 
(b) A = V2g being the initial state |+,0); (c) A = A\2g being the initial state |e,0); 
(d) A = being the initial state |+, v2i); (e) A = \/~2g being the initial state +, v2i); 
(f) A = 4\/2<? being the initial state |e, V2i). 



position of a nonrelativistic particle would do. 

While for the limits of zero and large mass, the state of the field remains Gaussian 
during the evolution, this is no longer true when the mass of the particle is such that 
mc 2 becomes comparable to c(P). This is the case of Figs. ITJ^b) and (e), where 
h\/2 = hg/y/2. Starting from initial states |+,0) and |+, \/2i), respectively, the 
Gaussian states of the field evolve into a nonclassical ones with negative Wigner function 
through a highly nonlinear interaction. The effect is more evident for the case of an 
initial state |+,0) in Fig. Hfb). Initial states with non-zero kinetic energy, such as 
|+,V / 2^) in Fig. (He) makes the qubit-field coupling term in Eq. flE]) dominant. This 
interesting finding shows that a Dirac evolution generates nonclassical features mainly 
in intermediate relativistic regimes, assuming that one traces out over the qubit degree 
of freedom. 

With our scheme we can also study the appearance of Zitterbewegung for a free 
spin— 1/2 Dirac particle. In Fig. [21 we plot the time evolution of (x) for a particle 
prepared in an initial state |+, 0) and with three different values of the mass. For better 
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comparison, plots (a) and (b) correspond to calculations done using the time dependent 
Hamiltonian ([6]) and the effective one (jSJ). The good agreement found indicates that 
choosing a value of the strong driving amplitude Q = 2n x 200 MHz is enough for our 
purpose. The dashed line in Fig. [2] shows the evolution of a massless particle, A = 0, 
whose expectation value of the position increases monotonically at a speed of light 
related to hg/\/2. If we assume a mass such that hX/2 — Kg/y/2, then we observe that 
the evolution of the (x) followed by dotted line departs considerably from the constant 
rectilinear behaviour one would expect for a free particle. Likewise, for a mass even 
larger, choosing for instance hX/2 = 4 x hg/y2, the effect seen in the dash-dotted line 
is even more significant. Now the particle hardly moves forward. Instead, it develops 
an almost stationary oscillation around its initial position. 

To gain in physical insight, we can derive an analytical expansion for the evolution 
operator of a free Dirac particle, 



At / /A 2 t 2 g 2 t 2 „ 

exp (—iHjyt/K) = exp I —ita y \j — + ^—p 2 J — io- z — sine \ — 1 — — p 2 



I g2j-2 




9 2 J\ 




-2 P ) 


— ia 




| sin | 


2 ) 




.At 




- 1— cr z smc( 



smc(gtp/ y/2) cos(gtp/ y/2) — sinc(gtp/ \/2) \ 

cos (gtp/ y/2) - sine (gtp/V2) , . 

- % — az ?tv + --" ( } 

where sinc(x) = sm(x)/x. Here, one can easily see that a larger initial average 
momentum gt(p(0))/2 will lead to a faster collapse of the trembling motion for a particle 
prepared in a wavepacket [HI [12] . 

As already discussed, Zitterbewegung stems from the fact that the time derivative 
of the position operator is not a constant of the motion for the Dirac equation. 
This can be seen if one writes down the solution for the evolution of this operator 
02]). The Zitterbewegung term Z(t) will be exactly zero if the initial state of the 
particle contains purely positive or negative components of the spinor. However, for a 
completely arbitrary initial wavepacket Z(t) does not necessarily vanish. In addition, the 
preparation of a purely positive spinor for a massive Dirac particle is not straightforward. 
In Ref. [15], a heuristic method for preparing a state with only negative components 
was used. This shows that even for a massive Dirac particle, evolutions without 
Zitterbewegung are possible. Another alternative to study the existence/absence of 
Zitterbewegung would be to implement a measurement of the Foldy-Wouthuysen position 
operator [9] |25] , 
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Figure 2. Expectation value of the a;— quadrature of the electromagnetic field mode 
inside the resonator. The interaction time is set to 30 nsec. The physical parameters 
used are g = 2vr x 10 MHz, Q = 2tt x 200 MHz, f = 0, together with: A = being the 
initial state |+,0) (dashed line); A = \/2g being the initial state |+,0) (dotted line); 
A = 4v2<? being the initial state |e, 0) (dash-dot); A = 4y2<? being the initial state 
|e, 0) using the FW Hamiltonian (solid line). Whereas in (a) the calculations involve 
time dependent Hamiltonians ([6]), in (b) effective ones (j8|) are used. 



X FW (t) = X + c 2 PH^H + A, (12) 
where A is the time independent operator 

s hca x hc 3 a x P 2 , , 

A = — ^ s — ^ , 13 

2E 2E 2 (E + mc 2 ) 



o~ x = i(?zO~y, E — v c 2 P 2 + m 2 c A . The Foldy-Wouthuysen position operator in (IT2|) 
has the property of being also canonically conjugate to the momentum operator. The 
main drawback of Xp^(t) from an experiential point of view is that is not diagonal in 
the coordinate space. This means that measuring such an operator might require a full 
quantum tomography of the evolved state of the system. 

Here, we follow a different route to the problem. Foldy and Wouthuysen are also 
credited for having introduced a unitary transformation [25] that allows for the direct 
diagonalization of the spin degree of freedom in the Dirac Hamiltonian. The physical 
relevance of this transformation is that it allows for a simple correspondence with 
classical mechanics [26J . Using the parameters of our simulation, this new Hamiltonian 
is given by 



h 2 \ 2 h 2 g 2 
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where 



S FW = exp j-i atari yj-j=a x pj j . (15) 

If one were able to follow the evolution of the system in the Foldy-Wouthuysen 
representation, but without transforming the measurement apparatus, then no 
Zitterbewegung would be seen. That is precisely, what we have done in the solid line of 
Fig. 12(b). This shows the evolution of (x) for a massive particle with hX/2 = 4 x hg/y/2, 
just like the dash-dotted one, but using the evolution operator 



Upwit) = exp (— iHwit/U) = Spw exp {—iHnt/K) Sp W (16) 

instead of U^it) = exp (—iH^t/h). Now the particle remains centered at (x) = 
during the evolution (solid line). More importantly, for large masses we can propose 
a physical implementation of the Foldy-Wouthuysen transformation using the same 
scheme we have developed. Clearly, when g(p)/(\\/2) is small enough, we can linearize 
the interaction in the unitary operator Sfw- This means that 



Z4w(<) ~ exp I -i-^=a x pj exp (-iK^t/h) exp \i^—^=a x p j , (17) 

where the first and last step involve the realization of a Dirac equation for a massless 
particle with an evolution time t = A -1 . Such an implementation would require some 
additional local rotations and Ramsey-like driving pulses |2l] . It would be particularly 
straightforward if the setup allows for a variable coupling from —g to +g, being g in 
strong-coupling regime. This is the case for instance, of a gradiometer type of coupling 
seen in flux qubits [27]. The solid line in Fig. |2(a) depicts the evolution of (x), computed 
under these assumptions, for the same particle of large mass, with h\/2 = 4 x hg/y/2. 
Even for such moderate value of the mass, it is evident how oscillations seen for the 
Dirac evolution become significantly reduced. 



4. Dirac particle in an external potential: Klein tunneling 

The addition of an external potential has the effect of coupling the different energy 
components of the spinor. As before, in the evolution under an external potential we 
can identify two limiting cases. A massless Dirac particle is described by the Hamiltonian 
Hk = hg/ y/2o~yp+h£\/2x. In Figs. [3] (a) and (d) the evolution of the initial state |+, 0) 
and |+, V2i), respectively, is shown. We use the same set of parameters as in Fig. CEJ 
with the only addition of a linear potential of strength £ = g/2. Again, starting with 
the qubit state, |+), allows for a simple interpretation. The initial state of the field 
remains coherent while being subject to two independent displacements. Namely, as for 
the free-particle case, there is one in the a;— quadrature proportional to g/y/2, and the 
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Figure 3. Wigner function W(x,p) representation of the state of the electromagnetic 
field mode inside the resonator. The initial state evolves under the Hamiltonian © for 
a time 60 nsec. The physical parameters used are g — 2ir x 10 MHz, ft — 2ir x 200 MHz, 
£ = g/2 and together with: (a) A = being the initial state |+, 0); (b) A = ^/2g being 
the initial state | + ,0); (c) A = 4\/2.g being the initial state |e,0); (d) A = being 
the initial state |+, V2i); (e) A = \[2g being the initial state |+, V2i); (f) A = 4\/2.9 
being the initial state |e, y/2i). 

second one induced by the potential that drags it down in the p— quadrature. Hence, 
the existence of such an external potential cannot alter the rectilinear movement in 
position representation at the speed of light. This is precisely the effect of the Klein 
paradox. A massless Dirac particle may propagate through the potential barrier with 
finite probability, regardless of the potential strength. The tunneling is accompanied by 
the 'rotation' of the positive energy components of the initial spinor. In our simple case, 
this means that the components with positive momentum components in the Gaussian 
wavepacket will be pushed downwards by the potential becoming negative. And the rate 
at which this happens is given by £ \[2 . This explains why the only noticeable difference 
between the plots in Figs. [Hand [3] is in the value of (pit)). 

If the mass of the particle is large enough, then the nonrelativistic Schrodinger 
equation should be a good approximation. Its Hamiltonian can be written as "HNRei = 
ha z p 2 /X + h£\/2x, which guarantees that any initial Gaussian state should remain 
Gaussian during the evolution. This is indeed what we see in Figs. [3][c) and (f) where 
we have prepared initial states |e, 0) and e, V2i), respectively, with a mass such that 
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h\/2 = 4 x %/v2- In the first of them, the particle is being pushed backwards by the 
potential. The same can be said for the second. Although having an initial positive 
kinetic energy has allowed it to enter the barrier further, after 60 nsec it is already 
moving backwards. 

While for these two limiting cases we see complete transmission or reflection, a 
particle with an intermediate mass will present only partial transmission/reflection. This 
is shown in Fig. [3]^b) which corresponds to the initial state |+,0) to a massive particle 
with hX/2 = Hg/\/2. The wave packet has now broken up into spinor components of 
different sign which move away from the center of the potential. If the wavepacket has 
some initial kinetic energy, as in Fig. E^e) where we prepare |+, y/2i), then the particle 
tunnels slightly to stop and break up eventually. 

As we did in Eq. ( fTTl) . one can obtain an analytical expansion for the evolution 
operator of a Dirac particle under the effect of a linear potential. Using the formula 



we see that for short time intervals the evolution corresponds to that of the free 
Dirac particle ( TTT1) . but with an initial state with rotated spin and displaced initial 
momentum with respect to the prepared wavepacket. 

The reconstruction of the Wigner function will allow to characterise what type of 
dynamics is being simulated. This can be done either using intracavity techniques [28] 
or allowing the quantum field to leak the superconducting resonator for its subsequent 
measurement. This will require the use of methods for reconstructing the moments of 
a propagating quantum microwave signal [2T| [22] . already available in cQED. 

5. Conclusions 

We have shown that circuit QED provides a powerful platform for simulating the 
relativistic quantum physics, establishing a useful dialogue between unconnected fields. 
Tuning the parameters of three classical microwave drives, the system is made to evolve 
according to the Dirac equation for a particle in presence of a linear potential. The 
degree of controllability offered by this scheme would allow to study interesting physical 
phenomena far beyond what any direct experiment has achieved. A relevant example is 
the possibility of implementing the Foldy- Wouthuysen transformation, which for decades 
has remained a purely abstract construction. 
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